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The effect of vertex corrections on the possibility of chiral symmetry breaking, 
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In this paper we consider the possibility of chiral (charge or spin density wave) symmetry breaking 
in graphene due to long-range Coulomb interaction by comparing the results of the Bethe-Salpeter 
and functional renormalization-group approaches. The former approach performs summation of 
ladder diagrams in the particle-hole channel, and reproduces the results of the Schwinger-Dyson 
approach for the critical interaction strength of the quantum phase transition. The renormalization- 
group approach combines the effect of different channels and allows to study the role of vertex 
corrections. The critical interaction strength, which is necessary to induce the symmetry breaking 
in the latter approach is found in the static approximation to be Oc = ^ /(cvf) ~ 1.05 without 
considering the Fermi velocity renormalization, and Oc = 3.7 with accounting the renormailzation 
of the Fermi velocity. The latter value is expected to be however reduced, when the dynamic 
screening effects are taken into account, yielding the critical interaction, which is comparable to the 
one in freely suspended graphene. We show that the vertex corrections are crucially important to 
obtain the mentioned values of critical interactions. 

PACS numbers: 


I. INTRODUCTION 

The fascinating physics of Dirac fermions, interacting 
via static Coulomb potential, became a subject of inten¬ 
sive research in condensed-matter physics after the ex¬ 
perimental realization of graphene^. Despite vanishing 
density of states at the Fermi level in clean graphene. 
Coulomb interaction may lead to different instabilities of 
electronic system^r.^. Even before the experimental dis¬ 
covery of graphene, the possibility of excitonic instability 
(related to the charge density wave formation) in this ma¬ 
terial was predicted theoretically^. It was argued, that 
sufficiently strong Coulomb interaction is able to open a 
gap in the excitation spectrum analogously to the chiral 
symmetry breaking in quantum electrodynamics^. 

The analysis within the Schwinger-Dyson mean-field 
approacb^rii in the static approximation found critical 
Coulomb interaction, which is necessary to open the 
charge gap, to be = e^/(euF) = 1.62 (see, e.g., dis¬ 
cussion in Ref. 0), where e is the critical charge, e cor¬ 
responds to substrate dielectric constant, and vp the 
Fermi velocity. This result is not very far from the result 
of quantum Monte-Carlo approacbi^, ac = 1.08. Inclu¬ 
sion of the dynamic screening of Coulomb interaction in 
Refs. [T^IT^ allowed to improve further the agreement 
with the numerical results. These results were however 
altered recently by considering the effect of renormaliza¬ 
tion of Fermi velocity, yielding critical Coulomb interac¬ 
tions, which depend on the approximation used, namely 
the critical couplings ac = 1.75 (Ref. [H), ac = 7.65 (Ref. 
[ 1 ), and ac = 3.1 (Ref. [ 13 ) were obtained. The latter 
result seems to be most reliable within the Schwinger- 
Dyson approach, since it accounts for the renormaliza¬ 
tion of quasiparticle weight and the polarization bubble. 
These later results are all larger than the coupling con¬ 
stant of free-suspended graphene, in agreement with the 


experimental fact, that freely suspended graphene shows 
no chiral instability^^. On the other hand, the important 
effect of the other bands for increasing critical Coulomb 
interaction was emphasized in Ref 

In the present paper we concentrate on studying the 
effect of the long-range part of Coulomb interaction 
and address the question, in which parameter range the 
Schwinger-Dyson approach is applicable. Indeed, this ap¬ 
proach, which can be reformulated as a ladder approxi¬ 
mation for electron interaction vertices, neglects contri¬ 
bution of the channels, other than the particle-hole one. 
From previous studies of the two-dimensional Hubbard 
model with short-range (on-site) interaction2i22ri2^, it is 
known that the contribution of the particle-hole channel 
is to a large extent compensated by the particle-particle 
channel. For the systems with Dirac electronic dispersion 
this cancellation is almost perfect, since the electronic 
Green function is odd with respect to momentum and 
frequency. At a first glance, this can strongly weaken the 
possibility of excitonic instability in graphene. 

Apart from that, the discussed approaches neglect 
electron-electric field vertex corrections, which also ap¬ 
pear beyond the ladder approximation, used in these 
studies. Although these corrections are expected to have 
no infrared divergencie a^'^i^^ , in the lowest (second) or¬ 
der in the Coulomb interaction they increase renormal¬ 
ized dielectric constant of the substrate. Therefore, one 
could expect that they yield further increase of the criti¬ 
cal value of the Coulomb interaction for the chiral quan¬ 
tum phase transition. 

The possibility of chiral phase transition in graphene 
was considered previously within the field-theoretical 
renormalization group approacb^^^— . It was shown that 
in the one-loop approximation^^— Coulomb interac¬ 
tion is irrelevant due to increase of the Fermi velocity 
with respect to its bare value. Apart from that, small 
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point-like interactions, which can be added to test the 
(in)stability towards chiral symmetry breaking, are also 
irrelevant because of the abovementioned cancellation be¬ 
tween particle-particle and particle-hole channels. There¬ 
fore, this approach does not yield divergence of some cou¬ 
pling constants at sufficiently large Coulomb interaction, 
and therefore does not yield chiral phase transition. In 
the two-loop approximation^ there is an instability of 
the field-theoretical RG flow at Oc = 0.78, which can 
be in principle related to a possibility of chiral phase 
transition. However, such a scenario is anyhow different 
from the two-particle instability considered within the 
Schwinger-Dyson approach. 

In the study paper we apply the functional 
renormalization-group (fRG) approacb^iSSvS^ and argue 
that considering momentum dependences of the vertices 
and summing all contributions from different channels in 
the one-loop approximation, the so far neglected effect 
of the vertex corrections yields at the same time enhanc¬ 
ing the tendency towards excitonic instability, and to a 
large extent compensates the above mentioned factors, 
which weaken the instability. Taken together, they yield 
even smaller result for the critical interaction strength, 
than in the ladder approach. In the static approximation 
without the renormalization of Fermi velocity we obtain 
ac = 1.05, which is smaller than previous estimates at the 
same level of approximation within the Schwinger-Dyson 
approach. The critical interaction increases to = 3.7 
when the Fermi velocity renormalization is taken into ac¬ 
count. The latter value of the critical interaction is how¬ 
ever expected to be reduced when the dynamic effects are 
taken into account, yielding the critical coupling, which is 
comparable to the one in the freely suspended graphene. 

The outline of the paper is the following. In Section II 
we formulate the model, in Sect. Ill we discuss reformu¬ 
lation of Dyson-Schwinger approach as summation of lad¬ 
der diagrams and extend previously obtained results to 
the symmetric phase. In Sect. IV we consider functional 
renormalization-group approach, discuss parametrization 
of the vertices used in the present study, and present main 
results for the electron interaction vertices. The Conclu¬ 
sion is presented in Sect. V. 


7 ^ are the Dirac matrices, e.g. 


7o = 



0 

-as 


7a = 




( 2 ) 


where ai are the Pauli matrices, a = Kevp) is the 
dimensionless coupling constant. The quadratic part 
of the action 0 represents the continuum limit of 
the microscopic tight-binding model of graphene (see, 
e.g., Refi^) and corresponds to the representation ip = 
'(/'bj V'B) “V'aI i'4’T is an annihilation operator for 
the electron in valley m and sublattice s). The scalar 
field Aq accounts for the Coulomb interaction of Dirac 
particles, which is integrated along the third coordinate 
z, transverse to the graphene layer. Excluding z direc¬ 
tion, perpendicular to the graphene layer, and performing 
Fourier transformation, yields the effective actioit^^^— 


La = - I Ipps [i7f,PtJ.^{q) + no^Oq] 'ipp+q,s 


(27r)6 

1 f (fq 
47ra ./ (27r) 




( 3 ) 


where = {v/vFiPxiPy) is the 3-vector, including Mat- 
subara frequency v (and similar for q). It is well known 
(see, e.g., Refi^), that due to vanishing density of states, 
the Coulomb interaction in graphene is not screened, 
since the particle-hole bubble 


n(q,a;) 


+ (w/'Cb)^ 


( 4 ) 


(where q = |q|) is itself linear in momentum at w = 0. 
Therefore, summation of the bubble contributions in a 
rundom phase approximation (RPA) yields renormaliza¬ 
tion of the dielectric constant only, such that the effective 
’’screened” Coulomb interaction is 


ar = a/{1 + 


( 5 ) 


III. THE LADDER (BETHE-SALPETER) 
APPROACH WITHOUT FERMI VELOCITY 
RENORMALIZATION 


II. THE MODEL 


To study electron interaction in graphene we consider 
the model (see, e.g.. Refs.— “— ) 

Lo = - j d^XiP^{-f^df, + i-foAo)iP, 

+ ^^^2 J d^^dz [(i9aAo)^ + {dzApPf''\ (1) 

where ip^ = ips{X) is the 4-spinor describing electronic 
degrees of freedom, s =t,4. is the spin projection, X = 
(T,x,y), p = 0,1,2, a = 1,2, a = {dr.VFdx.vpdy), ip = 


The Coulomb interaction may lead to excitonic (charge 
or spin density wave) instability with opening the gap in 
the electronic spectrum^iii2r2i. The simplest way to treat 
this possibility is the Bethe-Salpeter approach. Here we 
consider the instability from the symmetric phase, which 
complements previous studies of symmetry-broken phase 
within the Schwinger-Dyson approacb^i^rJi. The (lin¬ 
earized) gap equation of the Schwinger-Dyson approach 
in the static approximation ([5]) has the form 


Ap — Aq 


27r Ckr 


(M 


rTr 


7o ® 7o 


|p-k| 


nph(k, 0) 


ik, 


0 


( 6 ) 
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where we have introduced the bubble 


f /7p 

nph.pp(k, q) = - —[G{k,ie)^G{±'k + q,±ie)] 


= ± 


Vf 


2(uF(|k|)|k| +UF(|k±q|)|k±q|) 


ka{kb i Qb) 


( 7 ) 


where (p,k)> = max(p,k) and (p,k)< = min(p,k), k< = 
(k,k')<, k> = (k,k')>, Ki{x) = { 2 /tt)K{x), K{x) is the 
complete elliptic integral of the first kind. By considering 
the linear combinations 


'S 

kk' 

— kkl;/ -I- 2rj(i;, -1- F^^/ , 


\d 

kk' 

■pO OT~'^ _i_ 

— ^ kk' ~ kk' ^ kk' 5 


\a 

kk' 

pO pl2 

— ^ kk' ^ kk' 1 

(13) 


(the bubble Ilpp will be used later, in Sect. IV), G(k, e) = 
{jqE + mF(|k|) 7 ,jfca) , the matrices in the direct prod¬ 
uct act in the space of spinor indices of each fermionic 
line, Aq is the bare value of the gap, and in general the 
renormalized Fermi velocity z;F(|k|) depends on momen¬ 
tum. In the present Section we neglect the renormaliza¬ 
tion of the Fermi velocity by taking UF(|k|) = vp. The 
Eq. dH) can be iterated to express the resulting gap 


Ap = An < 1 


d^k 


Tr 


rpkonph(k, 0) 


( 8 ) 


through the vertex function rkk'q 7 for which we obtain 
the integral equation 


Stt cy^' 

-kk'q = n^_j^,| 7o^7o 

d^p 27raf, 


(9) 


+ 7o®7o 


(27r)2 |k-p 


■nph(Pj q)rpk'q) 


corresponding to the summation of ladder diagrams in 
the particle-hole channel. At q = 0 we therefore have 


kk'O — 


2710;, 

|k-k 


7i7o ® 7o + Y / / 


dOr, 


1 


27r |k — p| 


/®/-f ^^7o7a®7o7b ) Tpk'o- (10) 


we decouple these equations to obtain 


T^s,ci 

^ kk' 


r 


d 

kk' 


277 OLj' 

k> 


Ki 



^arrs,a f dp f J.s,a 

^ 4 y 



(14) 


with Ts = 2, Ta = 1. The solution to the resulting equa¬ 
tions can be searched in the form 


T^S,(1 

^ kk' 


k5«.a (k) , k>k', 
(^) 7 k < k'. 


(15) 


Often used approximation is to approximate Ki{x) by 
its small x asymptotic value, which is equal to 1. This 
corresponds to approximating 

l/|k-k'|«l/k> (16) 

in the equation (uni). The considering s-wave component 
for the vertex then coincides with the vertex itself, since 
the latter depends on the absolute values of momenta 
only. The solution to the Eqs. (HI reads (see Appendix 

A) 


We search the solution to this equation in the form 

Ekk'o = r^k'o7o®7o+rkk'o7a®7a 

+ rkk'o*7o7i72® *7 o7i72- (H) 


9m{x) 


\/l - UrTraX'^^ ’ 


(17) 


Performing matrix multiplication (see Appendix A), for 
the s-component of the vertices which are obtained 

by averaging over direction of vecotrs k and k' and de¬ 
noted (we use here a short notation k = |k| etc.), we 
obtain the system of integral equations 



7m = ^ (1 - VI - OirTm) ■ (18) 

Therefore, we find. 


p 0,12 

TTar 1 

1 Tra, 1 

TT Of.'p 

^ kk'O 

2V1 - 2a, k/k^-^ 

VI - ar k^^k^-^^ 

^ 2k>’ 

r),, — 

Tra, 1 

TT CXip 

(19) 

^ kk'O 

2V1 - 2a, k/k^-'’' 

2k>’ 

where 





7 = 7^,71= 7a- (20) 

The vertices diverge at a^c = 1/2 which coincides with 
the corresponding value, obtained from the analysis of 
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the symmetry-broken phase in Refiii^— . At 
we find 




CHr* 


p0,l,12 ^ 

^ kk'O — 


4V1 — 20 ^ (k<k>)i /2 


( 21 ) 


More accurate analysis, considering angular dependence 
of the Coulomb interaction, yields momentum depen¬ 
dences, which are qualitatively similar to the Eq. (HU), 
with the exponents 7 , 'Ji, determined by the Eq. (1451) in 
the Appendix A. The exponent 7 coincides in this case 
with the result of Ref. 0 , obtained in the fall on the cen¬ 
ter problem. The critical interaction in this case reduces 
to arc = 0.46, see Ref. E 

The chiral susceptibility can be obtained by perform¬ 
ing summation of the ladder diagrams for the trianfular 
vertex Ak via Eq. see Appendix A. To obtain the 
same result from the Eq. we will also need sublead¬ 
ing 1/A^ corrections to the vertex Tkk'o- Evaluation of 
these corrections can be performed in the way, similar to 
described above, and yields 


+ 2 rik'n + 


27Tar 


kk'O 


kk'O 


y/1 — 2oir 
7 


1 


k^-^k^ 


( 22 ) 


l-7Ai-27(kk')'' 


Integration of this vertex yields the triangular vertex (as¬ 
suming Aq = 1) 


Ak = 1 + i / (rL^o + 2rJ,k'o 


kk'Ol 


dk' 

27r 


2'y^l — 2ar 

1 A^ 

1 — 7 k7 ’ 


1 - 


AT' 


(1-7)2J kT 


(23) 


which is in agreement with the direct ladder summation 
(see Appendix A). We see that in agreement with the 
scaling analysis of Ref^ the triangular vertex is non¬ 
singular at the chiral phase transition; the singularity of 
the leading and subleading contributions to the vertex are 
cancelled in Ak. The corresponding chiral susceptibility 


X = 2 


dk 

— Ak = 

27r TT 


A 


A 47 
(1 — 7 )^ 27raf. 1 — 7 


(24) 


also does not diverge at the chiral phase transition. 


IV. FUNCTIONAL 

RENORMALIZATION-GROUP ANALYSIS 

A. General equations and parametrization of 
vertices 


The considered ladder approximation treats only a cer¬ 
tain set of ladder diagrams in the particle-hole channel. 


To consider the impact of other diagramatic contribu¬ 
tions, we apply in the present Section the functional 
renormalization-group approach^^i^i^^. This approach 
considers evolution of the electron interaction vertices 
due to different scattering processes. Since the contribu¬ 
tion of various channels of electron interaction can com¬ 
pensate each other (see discussion below), apriori the ap¬ 
plicability of the ladder approximation, considered in pre¬ 
vious Section, corresponding to the account of particle- 
hole diagrams only, is not clear. 

To treat the possibility of chiral symmetry breaking 
within the functional renormalization-group approach, 
we rewrite the action ® in the form 


Ln = -i 


+ I 




(27r)3 
d^pd^p'd^q 
(27r)9 


(25) 


, 27rQ! 


('0ps5'o'0p-t-g,s) I I ('^p^s^b'o'*/^p^ —g,s') 


and consider the Wick-ordered functional RG flow 
equatio n3^i32 for the electron interaction vertex 
Vi^ ^{ki,k 2 ,k^) (u, 2 ,fci ,2 and 13 , 4 , ^ 3.4 correspond 

to valley-sublattice indices and momenta/frequencies of 
incoming and outgoing particles). The corresponding 
equation can be written schematically as (see also Fig. 
1 ), 


14 = -Va o nA,pp o Va (26) 

-|-Va O IlA.ph O 14 

-I-2I4 * (nA,ph o —nA,ph*) Va, 


where (fc, q) = (k, (Tk -f q, + 

w) -I- F^^*^(k, j^)£))^*'‘(=Fk-I-q, -I- w) is the fermionic 
bubble, constructed from the Wick-ordering Green func¬ 
tions _^(k, v) and single-scale propagators J4(k, v) (cf. 
Ref. 1371), o and * stands for the summation over mo¬ 
menta, frequencies, and quantum numbers according to 
the diagrammatic rules (for explicit form of the flow equa¬ 
tion see Appendix B, Eq. ([5(1 1. 

The right-hand side of the Eq. (IMl) contains several 
terms, corresponding to the contribution of the particle- 
particle (pp), particle-hole direct (ph) and particle-hole 
crossed (phi, phi') processes (see also Fig. 1 and Ap¬ 
pendix B). It is important to note that due to the dis¬ 
persion, which is odd in momentum, the contribution of 
the particle-particle and particle-hole channels can com¬ 
pensate each other, since for small external momenta the 
corresponding bubbles in the first two diagrams in the 
right-hand side of Fig. 1 are equal in magnitude, but 
opposite in sign (see Eq. dT])); last three diagrams also 
yield some cancellation because of the sign change in the 
closed loop, and factor of two in the spin summation in 
the last diagram. 

These cancellations become transparent in the field- 
theoretical renormalization-group approacb^S^— , which 
was applied previously to the same model. In particular, 
in this approach the a^ contributions to renormalization 
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of short-range interactions, which would drive the chiral 
instability similarly to quantum chromodynamics^^^, 
are absent. At the same time, field-theoretical renor¬ 
malization group approach does not allow treating sin¬ 
gular momentum dependences, which are generated by 
the ladder diagrams in the particle-hole channel, if the 
abovementioned cancellation between particle-hole and 
particle-particle channel is weakly lifted by finite ex¬ 
ternal momenta. Therefore, we apply the functional 
renormalization-group approach to the considering prob¬ 
lem in this Section. 

In the pre sent study we use sharp momentum cutoff 
(cf. Refs.lm^) 


DA{k,e) 

FA(k,e) 


G(k,e) 


d(A-|k|) 
d(A-|k|) ’ 


(27) 


where G = (jqS + ivrJa^a — S(/c)) ^, neglect frequency 
dependence of the vertices, and either neglect self¬ 
energy effects (taking S(fc) = 0) or take the self¬ 
energy equal to its starting mean-field value, S(k) — 
(Q;/ 4 )uF 7 aA:Q ln(Auv/|k|), which coincides with the first- 
order perturbation theory result; Auv is an ultraviolet 
cutoff of the theory. The latter self-energy was recently 
shown to describe well the Fermi velocity renormaliza¬ 
tion in free suspended graphene, see, e. g., Refi^. The 
flow equations contain in both cases only the polariza¬ 
tion bubble, which is summed over fermionic Matsubara 
frequencies, 

nA,pp,ph(k, q) = ^ ^ l-lA,pp,ph(^; ^)aj —0 
iVn 

= npp,ph(k,q)d(|k|-A)0(A-|Tk + q|) 

+ npp,ph(Tk -f q, iq)(5(| ^ k -|- q| A) 

X 6 »(A-|k|) (28) 

where Hpp ph is given by the Eq. ([7]) with uji^dkl) = vp 
or UF(|k|) = uf( 1 + (a/4) ln(Auv/(|k|)). For E = 0 the 
functional renormalization group equations (1261) are anal¬ 
ogous to previously studied for the Hubbard model on the 
square- (Refs. [20l4^ and honeycomb (Ref. 0) lattices 
except that we use the continuum electronic dispersion 
and concentrate on the effect of the long-range Coulomb 
interaction. The use of the Wick-ordered approach allows 
us however to include easily the renormalization of the 
Fermi velocity; it also potentially allows to include in fu¬ 
ture effect impurities in this approach, see the discussion 
in Ref. 113. 

We parametrize the vertices 14 i,, 4 (ki, k 2 , ka) in terms 
of the average momenta kpp = (ki — k 2 )/ 2 , kph = 
(k 2 -I- k 3 )/ 2 , kphi = (ki -I- k 3 )/ 2 , and the momentum 
transfers qpp = ki-|-k 2 , qph = k 3 -k 2 , and qphi = ki-ks 
in the particle-particle and particle-hole channels. For 
convenience we decompose vertex functions into the con¬ 
tributions of the corresponding channels. Written explic¬ 



FIG. 1: Diagram representation of the flow equation (1261) . 
Circles denote vertices Va, solid lines without dash correspond 
to Wick-ordering propagators Da, while the lines with dash 
denote the single-scale propagators Fa- The lines inside ver¬ 
tices show direction of spin conservation. 


itly, our parametrization reads 

kii,,4,A(kl, k2, k3) = ^PP ~ Rphljqpp) 

+ ^^? 4 .A(kph + qphl, kph, qph) 

T V)i ^A(^phl,kphl qph,qphl). 

(29) 


The first two arguments of functions correspond 

to the average incoming and outgoing momentum, while 
the third argument denotes the momentum transfer. 

The functions three continuum 2-component 

variables is hard to treat accurately numerically. For phi 
channel we pick out renormalized Coulomb interaction 


^i.. 4,A^^P^1’Qph,qphl) 

^i.. 4,A^^P^1’^P^l ~ Qph,qphl) 

27ra 


(30) 


(qphl) .Uplil ^A.kphi qphiUphl 


by introducing the renormalized polarization n(q) and 
electron-electric field vertex 5 A,k,q functions. We then 
approximate the functions ^ii^4 A neglect¬ 

ing their dependence on their third argument (i.e. 
transfer momenta qpp, qph, or qphi), which are as¬ 
sumed to be zero in the actual calculation of these 
quantities. The justification for this approximation 
is that at finite k, k' the dependence of these func¬ 
tions on qpp,ph is non-singular, while there is an es¬ 
sential singularity of the considering functions on the 
first two momenta, ^(k, k', 0 )^ max(A, |k|, |k/)^~^ 
X max(A, min(|k|, |k/))~^. We have also verified that 
taking into account the full momentum dependence of 
^i^4^a ^1^4 A yields only small corrections to the 

obtained results. 

To treat accurately the remaining momentum depen¬ 
dences, we follow the idea of Ref,^, expanding these 
functions in some basis. Since the dependence on the 
absolute value of first two arguments is expected to be 
singular (as follows from the analysis of Sect. II), we ex¬ 
pand in Fourier harmonics with respect to the angle of 
each of the two momenta. The resulting flow equations 
for b/PP^P^’ ^ii*'4.A> 5 A.k,q, and nA(q) are presented in 
the Appendix B. 
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A 


k 


k 

k 


FIG. 2: Diagram representation of the flow equation for the 
triangular vertex Ak.A (denoted as triangles in the right-hand 
side). The first (both) term(s) in the right-hand side are re¬ 
spectively included for the spin (charge) vertex. The nota¬ 
tions are the same, as in Fig. 1. 




To calculate the chiral spin or charge susceptibility, 
we first calculate the flow of the triangular vertex Ak.A 
according to the diagrams of Fig. O The susceptibility 
is then obtained straightforwardly by convolution of two 
triangular vertices with two Green functions. 

The flow is started at A = Auv > 1 /uf, and finished 
for some A = Amin ^ vf- The bare value of the vertex 
corresponds to the bare Coulomb interaction 

ni2 

aAu^.k.q 

nA„„ (q) 

GPP , 

^^1..4,Auv 

A il A2 
^k,Au^ 

B. Results 


(70)0*2 j 

|qL 


= 


U..4,Au 


= K' 


phi 


U..4,Au 


= 0 , 


= <5,. 




(31) 



FIG. 3: Scale dependence of the absolute value of the ver¬ 
tex functions |k||k|) (a,c), the renormalized 
Coulomb interaction aij,A(q), and the electron-electric field 
interaction vertex ^ (b,d) for a = 4/(27r) (a,b) and 
a = 5.5/(27r) (c,d). Vertices pp,ph,phl correspond to the 
intrasublattice interaction ii ,.4 = 1, while pp',ph' correspond 
to the intersublattice interaction ii = is = 1, is = ii = 2. 


1. Flow without the Fermi velocity renormalization 


First we neglect the self-energy effects by putting E = 
0. In the absence of self-energy effects, all the vertices 
have dimension k~^, and can be expressed in terms of 
the scaling functions 


nA(q) = 


^ rPP,ph,ph' / 1^ I ^ -^uv 

(32) 


I I ,n / A Aa 
'iql’lql 


i^A.k.q 


= f, 


2i ,i2 ,m,n 


Iql’lql’ |q| 


where to, n numerates Fourier harmonics with respect to 
the directions of momenta k, k', or q. For sufficiently 
large Auv I/'Cf and small momenta, one can ignore 
the last argument in these scaling functions. 

The flow of GP^(k,k') alone (with zero interactions 
\/pp.phi ^ reproduces the results of the ladder approach 
of Sect. II. We obtain the critical value a);’' = 3.0/(27r) 
which is close to the value = 2.9/(27r), obtained in 
Refsi^Aii, the difference is related to the discretization of 
momentum dependence of the vertices. The flow of n(q) 
alone reproduces (in the end of the flow) the renormaliza¬ 
tion of the dielectric constant by the static polarization 
bubble, described by the Eq. ®, 


f^{0,oo) = ZD = l + TTa/2. (33) 


In the present paper we consider however the com¬ 
bined contribution of all the channels. To analyze the 

results, we plot dimensionless vertices = 

|k|I7PP’PGP^^(|k|, |k|) as functions of A/|k|, coupling con¬ 
stants afl;,A(|q|) = a|q|/nA(|q|), and the vertex func¬ 
tions 5 a oq as functions of A/|q| (see Fig. |31). We see, 
that the scaling forms (15^ are approximately fulfilled. 
For small a < 5/(27r) we find that the interactions, gen¬ 
erated in the particle-particle and particle-hole channel 
are close to each other in the absolute value, having the 
opposite signs, and therefore, almost cancel each other. 
At the same time, for larger interactions, this cancellation 
is lifted, and one of the coupling constants, correspond¬ 
ing to the intrasublattice interaction in the particle-hole 
channel, becomes bigger than the other interactions. 

The interaction dependence of the in¬ 
verse dimensionless quantities Gq m n = 

lima;^o.y^oo(l,a;,y), corresponding to the 
limit A:—>'0,A/A:—^-Oof the coupling constants of Fig. |31 
is shown in Fig. 01 We find that similarly to the ladder 
approach of Sect. Ill, the leading intrasublattice compo¬ 
nent GoAiii (which is equal to |k|(r°,^(,-bri^o) — 2Trar = 
7 rar(I/-\/I — 2Q;r — I) in the ladder case), increases faster 
than other components. By fitting the obtained re¬ 
sults with the dependence Ayjyja^y —« we find 
OLc,v = 6.9/(27r). From this fit we also observe, that 
the ladder approximation with renormalized parameters 
is applicable in fact only at a > 5/(27r), i.e. outside 
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FIG. 4: (Color online) Interaction dependence of the absolute 
value of inverse small momentum s-wave components of the 
dimensionless vertex functions a-nd inverse charge 

susceptibility in the end of the flow. Long-dashed line cor¬ 
responds to the fit of ph vertex function by the dependence 
Av!y/cxc,v — a, long-dashed line corresponds to the ht of the 
inverse susceptibility by x~^ = ^x(l + The 

other notations are the same as in Fig. [S] 



FIG. 5: (Color online) Interaction dependence of the renor¬ 
malized Coulomb interaction q:h,o, the electron-electric field 
interaction vertex (a), and the exponents of the momen¬ 
tum dependence of the vertices Xvx shown by points. 

Upper solid line in (a) corresponds to the weak-coupling re¬ 
sult m of Ref.—, while lower solid line corresponds to the 
linear dependence 1-1-a/2; solid line in (b) is the result of the 
fit coording to the Eq. (I35| 


the region, where particle-particle and particle-hole 
channels compensate each other. The corresponding 
region below the critical interaction, which can be 
associated with the critical regime, appears to be rather 
narrow. Calculations of the chiral susceptibility shows 
that the spin and charge susceptibility are equal in 
the considering model, since the second contribution 
in the right-hand side of Fig. vanishes. Similarly to 
previous studies^i^, additional short-range interactions 
are expected to remove this degeneracy, yielding either 
spin or charge order depending on the ratio of the on-site 


and nearest-neighbour Coulomb interaction. Since the 
nearest-neighbour Coulomb interaction is expected to 
be smaller than the on-site component, the spin-density 
wave is expected to be more favourable, than the charge 
density wave. The susceptibility, obtained within the 
present model, containing only long-range interaction, is 
not singular near the chiral phase transition, similarly to 
the ladder approximation. Its fit with the dependence 
X~^ = + ^xV^c,x ~ similar to the obtained 

in Eq. (1^ . yields ~ 6.7/(27r), which is somewhat 
smaller Ucy- The small difference between Ucy and 
ac,x likely occurs because of the narrowness of the 
region, where such fits are applicable. 

We see that the critical coupling constants acy,x 
approximately twice larger, than the constant arc, ob¬ 
tained in the ladder analysis. To understand the reason 
of this difference, we plot in Fig. [5^ the renormalized 
coupling constant where 

Z'j^= lim nA(|q|)/|q| = hm f^{x,y), 

q^0,A/q^0 a:->.0,y->oo 

and the vertex 

= lim (z, X, y ). 

x,z—¥\j,y—¥oo ’ ’ ’ 

Due to vertex corrections, we obtain Z/, < Zr, which 
agrees with the observation of Refi^i^i that the vertex 
corrections increase the effective dielectric constant. The 
obtained behavior of Z'j^ agrees well with the result of 
Ref.^, 

Z/j = Zz,-b 0.78^2-bO(a3), (34) 

although the latter is obtained in the weak-coupling limit 
(second order in a). We also find, that the interaction 
dependence of the intrasublattice component is well 
fitted by g}^ ~ 1 -|- a/2. This linear dependence of the 
vertex (with the correct coefficient) can be obtained by 
comparing linear and quadratic terms in Eq. (IMl) . 

From Fig. [5] we find that near chiral quantum phase 
transition the coupling constant renormalization factor 
<ARfi/a ~ 0.3, which is smaller than the ratio of the 
ladder and fRG critical couplings arc/o-c — 0.43. This 
difference can be attributed to the vertex corrections. 
Although some vertex corrections (yielding suppression 
of the Coulomb interaction) are already accounted by 
Eq. (IMl) . one should take into account, that when con¬ 
structing a particle-hole ladder, in the presence of the 
vertex corrections each Coulomb line acquires a factor 
g^. We find, that (gj^)^ ~ 2.25 near the chiral phase 
transition. This itself would yield two times decrease of 
the critical interaction, which would be therefore equal 
Q;rc/0.3/2.25 ~ 0.68, somewhat smaller than the ob¬ 
tained value ac- The remaining difference can be ex¬ 
plained by partial compensation between particle-hole 
and the particle-particle channel near quantum phase 
transition. 

From Fig. [3J: one can see that at a close to chiral 
phase transition, the dimensionless coupling constants 











do not saturate yet in the considered range of A/q (the 
same behavior is observed for susceptibilities, which are 
not shown). Considering smaller A requires increasing 
number of discretization points, and not feasible numeri¬ 
cally. Apart from that, some deviations from scaling are 
observed, which may be attributed to the discretization 
procedure. To determine critical Coulomb interaction 
more accurately, we determine the exponents 7 ^ and 7 ^, 
defined by 

r.A,o.o,»(|k|,|k'|) = k'<K<k, 

Ak.A = k<k (35) 

A 'A 

These exponents can be determined in the whole range 
A <C 1 and do not require achieving saturation of the cou¬ 
pling constants. In the ladder approach one has 1 — 7 ^ = 
7 ^ = 7 . The results of the renormalization-group cal¬ 
culation are shown in Fig. [5}). We see that the ex¬ 
ponents 7 y, 7 ^ also behave similarly to the ladder ap¬ 
proach for 2TTa > 5, and show a bifurcation point at 
cxc = 1.048 ± 0.008. We consider this result as an es¬ 
timate for the critical coupling strength. This is much 
smaller than the value Uc = 1.62 in the ladder approxi¬ 
mation, and agrees well with the results of Monte-Carlo 
analysisi^. 

2. Flow with the Fermi velocity renormalization 

To model the effect of Fermi velocity renormalization, 
we put T,{k) = (a/ 4 )tiF 7 afca ln(A„„/|k|). The scale de¬ 
pendences of the vertex functions are shown in Fig. [51 As 
one can expect, the universality of the scaling functions is 
violated by renormalization of the Fermi velocity. Due to 
relatively large Coulomb interaction, the contributions of 
the particle-particle and particle-hole channels are split, 
such that the above discussed compensation of the two 
contributions is less pronouced. The electron-electric 
held vertex renormalization achieves values 7 ~ 1.5, 
which are comparable to those obtained in the previ¬ 
ous subsection without Fermi velocity renormalization 
(although for larger interaction strength), strongly sup¬ 
pressing the critical interactions. The dependence of the 
interaction in the end of the how on the coupling con¬ 
stants is shown in Fig. [T] We obtain the critical coupling 
constant a” = 3.7, while without the vertex corrections 
in the Bethe-Salpeter approach the chiral instability does 
not occur in the static approximation with the Fermi ve¬ 
locity renormalizatioiii^. Apart from the vertex correc¬ 
tions, discussed in previous Section, the reason of getting 
hnite and not large critical interaction is in the renormal¬ 
ization of Fermi velocity also when calculating the polar¬ 
ization bubble. This effect was not accounted in Ref. [l^ 
and reduces the screening of the interaction, decreasing 
therefore the ehective interaction. To see more explicitly 
the ehect of changing screening in the presence of the 
renormalization of the Fermi velocity, we switch off the 



FIG. 6: Scale dependence of the absolute value of the ver¬ 
tex functions |k||k|) (a,c), the renormalized 
Coulomb interaction a_R,A(q), and the electron-electric field 
interaction vertex (?a 0 q (b,d) for a = 16/(27r) (a,b) and 
a = 20/(27r) (c,d) with account of the Fermi velocity renor¬ 
malization. The notations are the same, as in Fig. [S] 



27ia 


FIG. 7: Interaction dependence of the absolute value of in¬ 
verse small momentum s-wave components of the dimension¬ 
less vertex functions Foqi, 1 , 1 , 1 }, 0,0 (line and points) and in¬ 
verse charge susceptibility (crosses) in the end of the flow with 
account of the renormalization of the Fermi velocity. Long- 
dashed line corresponds to the extrapolation of interaction 
dependence of the inverse vertex. 

renormalization of the Fermi velocity in the polarization 
bubble nA(q) and obtain = 4.6 which is somewhat 
higher than the a", but yet finite and even considerably 
smaller than the result of Ref. [l^ with dynamic effects 
included. Therefore, we conclude that even with the 
Fermi velocity renormalization, vertex corrections yield 
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substantial decrease of the critical interaction. The ob¬ 
tained value a" = 3.7 is expected to be further reduced 
by the effects of dynamic screening of Coulomb interac¬ 
tion, which may yield the critical Coulomb interactrion 
comparable to the one in suspended graphene. 


V. CONCLUSION 


In the present paper we have analyzed the results of 
the Bethe-Salpeter (ladder) approximation and the func¬ 
tional renormalization-group approach, which accounts 
for all the channels of electron interaction to describe 
chiral phase transition in a system of Dirac electrons. 
We have shown, that without the Fermi velocity renor¬ 
malization, at sufficiently small coupling constants a < 
5/(27r) the particle-particle and particle-hole channels 
partly compensate each other. On the other hand, for 
a > 5/(27r) the compensation is not present, and the 
particle-hole channel dominates. This yields chiral phase 
transition at ac = 1.05, which properties are rather sim¬ 
ilar to those, obtained in the Schwinger-Dyson or Bethe- 
Salpeter (ladder) analysis. The vertex corrections en¬ 
hance the tendency towards chiral symmetry breaking 
and compensate the effect of partial cancellation between 
the particle-particle and particle-hole channels. 

With accounting the Fermi velocity renormalization, 
the abovementioned compensation is weakened, and 
therefore vertex corrections are expected to have even 
stronger effect on the critical interaction sterngth. In par¬ 
ticular, in the static approximation we obtain a)! = 3.7, 
while this value is expected to be further reduced by 
the effects of dynamic screening of Coulomb interaction, 
yielding the results for the critical Coulomb interaction 
smaller than previous estimates. Whether it remains 
above or below the experimental a = 2.2 and how much it 
differs from the results of Monte-Carlo analysis, requires 
additional studies. As also described in the introduction, 
the effect of the other bands, which add short-range in¬ 
teractions on the top of the long-range Coulomb tail, can 
be also important. Their analysis is beyond of the scope 
of the present work, and is interesting to be performed 
in the future. Finally, the effect of impurities, either 
neutral, or charged, in the presence of the interelectron 
Coulomb interaction, has to be analyzed. 
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Appendix. A. Solution of the Bethe-Salpeter 
equation for the vertex 


Substituting the ansatz (HU into the Eq. (nni) we find 


X 


(rkk'o - 70 ® 7o + rkk'o7a o 7a (36) 

+rkk'o*7o7i72 O *7 o7o7i72 


Ck.j' 

T 


dpdOp 1 
27r Ip — k| 


7 0/ 


P<^Pb ^ 
-|^7o7a O 7o7a 


[rpk'o7o0 7o+rpk'o76 0 76 
+ rpim*7o7i72 O *7 o7i72] ■ 


In the following we consider the ”s-wave” component of 
the vertex, which is obtained by averaging over the di¬ 
rections of k and k'. For this component, which we refer 
in the following as F^,, we obtain. 


n 27ra; 

- -7—^1 

k> 


+rfcfc/*7o7i72 ' 


T 


dp 


{P, k): 


-Ki 


t< 

k> 

*7o7o7i72 

{P,k)< 

{P,k)> 


7o®7o + rL'7a®7a 


(37) 


I - 27o7a 


>7o7o 


X [rpfc/7o O 7o + rpfc'76 ® 7b + rpfc/l7o7i72 ® i7o7i72] , 


where Ki{x) = {2/Tr)K{x), K{x) is the complete ellip¬ 
tic integral of the first kind, fc< = min(fc, fc'), fc> = 
max(fc, fc'), {p,k)^ = max(p, fc), (p,k)^ = min(p, fc). 
Performing matrix multiplication, we obtain Eq. dni 
of the main text. 

With the anzatz dTSl) for the linear combinations (1131) . 
we obtain for the functions gm{x) with in = s, a the 
equations 




Oij.Tm 

Akk' 




(38) 


Supposing the last integral is convergent at A —>■ oo, these 
equations can be rewritten as 


gm.{x) = 27rar7fi {x) + J dpgm{p)Ki (px) 

dp 


OL'p'^m 

OirTm 

Ax 


— gmij^Ki 
p \p 


dpg 

m ip)Ki 


(39) 


The obtained integral equation can be then transformed 
to the differential one. 


x[xgm{x)]” 




gm.{x) 


(40) 


4 
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with 


gm{l) - 2Tra = -2g^(l) = dxgm{x). (41) 

The solution to the obtained equations has the form of 
the Eq. (I17L of the main text, which yields 


^>^kk’ — 


k r°d 2 _ 

K>i 1 . 1.1 — 


TT CXr 


ky\'^ irar 


2V1 - 2ar \k<J 2 ’ 


(42) 


TT CX'T 




^ 2V1 - 2a, V^<. 

k~^ \nar 


± 


TT Oir 


\/l — Olr \ fc< 


where we have denoted 7 = 7s = 5(1 — v^l — 2ar) ■ 

The solution of the equations (ESI) beyond the ap¬ 
proximation (HU can be obtained analytically only for 
X = k'/k <C 1. Assuming again gm{x ^ 1) = AmX~^^ 
we find: 


the result for the exponent 7 = 7 s coincides with that 
obtained in the fall on the center problentiS. For the 
constant Am we obtain the equation 


^r' 


arTr, 


47m 

ar'Tr) 


dpgmip) 


/o p L 


A 

/—A -^m 


= 27rQ; (46) 


which can be solved numerically. Qualitatively, the solu¬ 
tion do not differ from the simplified result (l42l) . 

Similar consideration can be performed for the trian¬ 
gular vertex, considered in Refj^, cf. also main text. 
The corresponding equation for the ” s-wave” component 
of the vertex reads 


Ap — Ao -I- — 


dk I 


(kiP) 


{k,p)-: 


A„iX = 27rQ; + 


axr. 


4 Jo 

arm dp 


dpgmip) 


(43) 


4x 


—gmip)Ki 
P \P 

[ d^~^^ Ki 


To simplify the second integral, we introduce the con¬ 
stant C such that a; <C C <C 1. Splitting the regions of 
integration over p by the constant C, we find 


AmX = 27rar -I- 


arTr. 


arX m,A. 


4 JO 

rC 


[ dpgmip (44) 

Jo 


dp / X 

P 


^+1 


OirTr] 


dp 

^ Jc P 

k^rXmAm f 
4a; ./n 


Pm ip) 

dpp~'^’^ Ki 


where Aq is the bare value of the triangular vertex. 
Substituting Ap = Ap~^ for p <C A and introducing 
p <C C -C A, we obtain for 7 the equation 


, -V A cirA f dk f k 
0 


(48) 


arA 


c 


A 


f dk 

J\ Oir [ 

1 

-^yJ 


— An 


P 

arA 

2pT' 

A 


dy 


c 

00 


Ki iy) 


dk 

dk 


Lo 


y 


1+7 


Ki 


ar f dk / A \ arAPA 

2 j k \ ^ k"<) 27 

0 


From this we find again 7 = 7 s, determined above, and 
the constant A is determined by 


arVm f i~ /—, OirPmAm 

= 2Trar -\ -^— / dpgmip) - 


OLr^r) 


47r, 


-b 


_ ^ ^ 

4 Jo p : 

^rXmAm 


A 

/—A -^m 


f dyy iy) 

Jo 


+ 


h y 


dy ^^.71 

■ y 


1 - 1-7 


The resulting equation for jm reads 


1 = 


■ f dy{y^-- ^ +y ^™) AT (p); 
^0 


(45) 


arAA^ ar 
2 


27 


dkf A 

T ^'■-F 


= An. 


(49) 


With the approximation ((T5)) this reduces to the result 
of the main text. 


Appendix. B. Functional renormalization group 
equations 

In this Appendix we present analytical form of the 
renormalization-group equations, shown in Fig. 1 of the 
paper. The equation for the vertices Td™^(ki, k 2 , ka) in 












































11 


the static approximation reads: 

T^n..4.A(ki,k2,k3) = - 

k — CJpjj, k)14^i2Z3i2,A(k, k2, k3) 

“t” ’ k) 1 ^^ 2443 44,A(k, k + Qpp, k3) 

xnl^;r*kk,qpp)+ 

k - Qphi, k 2 , k - Qp^) 

“t“ ^^12^232^ ,A (ki, k — Qphl ; k 3 ) 1^^222422 ,A (k, k2, k2 ~\~ Qphl) 
23 2^^A(ki, k qphi, k 3 )V^^i 2 22 ^ 4 ,A(^i ^27 k Qph)] 

xnX^*3k,qphi)}. (50) 

where kpp = (ki — k 2)/2 and kph = (k 2 + k 3)/2 are the 
average momenta in the particle-particle and the particle- 
hole channels, qphi = ki—ks and the momentum transfer 
in the latter two channels, qpp = ki - 1 - k 2 and qph = 
k 3 —k 2 . Using decomposition (IMI), we obtain the equation 
for the contribution of the particle-partice channel 


U,PP 3 k,k',q) = 

n..4.k" 


(51) 


^4^4i4'(k,k",q) 


u.p’" 

2 i 22 2 -|^ 22 


+ K 


phi 


21222 j ^22 


k-k" + q k"-k + q 

2 ’ 2 
k -f k" + q k + k" - q 


,k-|-k 

k-k 


U.PP. . (k",k',q) 

2 o 2 4 2.q24 '' ‘ ’ 


+ U?') 


k"-k' + q k'-k"-fq 


23242324 ^ 2 

/ , 1 '' 


+u.?v. . 

23242324 


k'-fk" -F q k' -f k" - q 
2 ’ 2 ' 

and the particle-hole channel, 
iAP‘' 3 k,k',q)= ^ 

7 ' Ic" 

k-k" + q k"-k + q 


k" - k' 


(52) 


u.pp,. 

'^1 ^ 4^1 ^4 
rph 


,k-f k' 


+ r"“;.;„(k,k",q) 

+v-p“ ,!S±!L±3,li±hU3,k_k'' 


41144444 


U-PP . 

43434342 


2 2 
k"-k'-t-q k'-k" + q 


+ r.F“.. 4 (k".k',q) 


+yp‘’P 


23 22 ^3 ^2 


k'-fk"-I-q k'-fk"-q 


k"-k 


We represent the vertex according to the Eq. (l30l) 
of the main text, and in the most part of calculations, 
we neglect the dependence of UPP’P*' and Up^^ on the 
third argument (as discussed in the main text), project¬ 
ing it to zero. We have verified that treatment of the 
full momentum dependence does not change the results 
substantially. For the vertex we obtain the equation 


^phl _ ^ 

Jk,k',q)= ^ nl=^*^(k",q) (53) 

*'i..4.k" 

/ k-k"-fq k-k"-q 

4i4i4343 2 ’ 2 ’ ^ 

+ V.P; , j k + V' + q k + k"-q 

212^2323 ' ^ t 


+h^E4'(k,k",q) 


/ k"-k' + q k'-k" + q 

4^42444^ y 2 ’ 2 ’ ^ 

^ph 


+ h?E4' (k",k',q) 

phi ,k" + k' + q k' + k"-q 


-, T/pni / 

^ ^ 44 * 24442 ' 


2 ’ 2 
k-k"-fq k"-k + q 


,k"-k': 
k-f k' 


U.PP,. 

2 l 2 i 23 23 

+ r.fU(k,k”,q) 

^phi / k-f k" + q k" -f k - q , _ , „ 

+ *^ 212 ^ 2^23 2 ’ 2 ’ 

UPP, f k"-k' + q k"-k'-q ^^ , 

24222224 \ 2 ’ 2 ’ 

+ UP^,. (k"+k' + q k" + k'-q 

44424^44 


+UP“,. (k",k',q) 

24222^24 ^ 


- 2 


| k-k" + q k-k"-q 

414x4343 (22 


H ('k + k" + q k + k"-q _ 

212^2323 ' ^ 7^7 


+hS34^(l^A",q) 

/ k"-k'-fq k"-k'-q 

'' 4 ^ 424^44 I 2 ’ 2 ’ ^ 


ph I k"+k' -f q k" + k' - q 


-f vr ■ 

24222224 


k" - k' 


+UP“,. (k",k',q)j|. 

*4424244 ' ’ 
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Finally, the vertex 7 and the inverse propagator 11 are 
determined from 


9T 


’X= E n;=X^nk",q) 

j' V" 

yPP /^k-k" + q k"-k + q 

il 2 ^ ig 23 


(54) 


,k + k' 


numerically. To treat accurately these dependences, we 
follow the idea of Ref.—, expanding these functions in 
some harmonics. Since the dependence on the absolute 
value of first two arguments is expected to be singular 
(as follows from the singular behavior of the gap func¬ 
tion A(k) ^ |kp^, 7 ^ 1/2 on the ordered side of the 
transition, cf. Refs— dd), we expand in Fourier harmon¬ 
ics with respect to the angle of each of the two momenta: 




-f 


2l 2-|^ 2g 23 


- 2 


y.pp. 

* 1 ^3 ^3 


k + k" + q k" + k-q 
2 ’ 2 ’ 
k — k” -f q k k” — q 


RP*) . 

2l 2-|^ 23 2g 


2 ’ 2 
k-t-k" + q k + k"-q 


k-k 

,k-f k" 


k-k' 




and 


n(q) = -2 ^ ngf^'^(k'',q)4;;^g^4. (55) 


The equations (IST]) - (I55|) have to be solved numerically. 

Yet, the functions and « of two contin- 

uum 2 -coinponent variables is hard to treat accurately 


np 

C™.4.A(k,k') = ^ [f//™’'^( fc,/c')cos(TO(^k)cos(m'(/3k') 

m—0 

+ G™™ k') cos(mv5k) sin(m'(pk') 

+ k') sin(mv9k) cos{m'ip^,) 

’^{k, k') sin(m:^k) sin)™':/?!,') 

(56) 

and discretizing the absolute values of the momenta, with 
further performing (bi-)linear interpolation between the 
discretization points. In this way, we achieve sufficiently 
fine discretization of each vertex, (typically we take np = 
2 Fourier components, corresponding to 2np -1-1 = 5 
Fourier harmonics, and n = 20 25 radial points log¬ 

arithmically distributed in the range [Amin/lO, Auv]; we 
take Auv = 2 and Amin = e“^°). The results of the solu¬ 
tion of Eqs. (l5T|l - (|56|l are discussed in the main text. 
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